%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% variance decomposition for structural identification
% this version: October 2010
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [frac_rev_med,frac_rev_low,frac_rev_high] = ...
    vardecomp(Respext,num_shocks,num_vars,HORIZON_ident,ndraws)

var_rev = zeros(num_vars,num_shocks,HORIZON_ident,ndraws);
sum_rev = zeros(num_vars,HORIZON_ident,ndraws);
frac_rev = zeros(num_vars,num_shocks,HORIZON_ident,ndraws);
for n = 1:ndraws
    % iteration
    dnum=num2str(n);
    if dnum(end)==num2str(0) && dnum(end-1)==num2str(0)
        disp('iteration'),disp(n)
    end
    for time_counter = 1 : HORIZON_ident
        for k = 1:num_vars
            for j = 1:num_shocks
                var_rev(k,j,time_counter,n) = sum(Respext((j-1)*num_vars+k,1:time_counter,n).^2);
            end;
            sum_rev(k,time_counter,n) = sum(var_rev(k,:,time_counter,n));
            for i = 1:num_shocks
                frac_rev(k,i,time_counter,n) = var_rev(k,i,time_counter,n)/sum_rev(k,time_counter,n);
            end;
        end;
    end;
end;
% sorting; note that this induces that the fractions do not add up to 1
frac_rev_sort = zeros(num_vars,num_shocks,HORIZON_ident,ndraws);
for k = 1:num_vars
    for time_counter = 1:HORIZON_ident
        for j = 1:num_shocks
            frac_rev_sort(k,j,time_counter,:) = sort(frac_rev(k,j,time_counter,:),4);
        end
    end
end;
frac_rev_med = frac_rev_sort(:,:,:,fix(0.5*ndraws));
frac_rev_low = frac_rev_sort(:,:,:,fix(0.16*ndraws));
frac_rev_high = frac_rev_sort(:,:,:,fix(0.84*ndraws));
